Jason 2 Two Dimensional Fourier Analysis

import sys
sys.path.insert(0,'..')
import Jason2_py as Jason2
import numpy as np
import pandas as pd

# Plotting Functions
import holoviews as hv
import cartopy.crs as ccrs
import plot_funcs as pf
from IPython.core.display import HTML

Plotting Modules:
GeoViews 1.6.6
HoloViews 1.13.3

Jason2_py

The Python package written for this project, Jason2_py, can be found at https://github.com/theo-yang/jason2. The package contains the following modules:

• Jason2_py.download.py: tools for downloading and reading Jason-2 ssha data directly from the online repository, as well as functions for reading netcdf files from Copernicus datasets. ECCOv4 datasets can be read using ecco_v4_py

• Jason2_py.preprocess.py: tools for removing ssha datasets based on NaN and distance thresholds.

• find_2D_power_spectrum.py: function for calculating the Hann-windowed two dimensional Fourier transform of ssha data.

• Jason2_py.model.py: tools for modeling the two-dimensional Fourier transform of ssha, as discussed here

See full dosctrings using ? operator

Datasets

The following analysis uses datasets from the following sources:

• Jason-2 for sea-surface height anomaly
• ECCOv4 for stratification climatology
• Copernicus for sea surface currents

Now let's download the data for this demonstration (you can also just download from saved data files):

Jason-2

# get existing cycles:
f = open("cycles.txt", "r")
cycles = [line.strip()[1:4] for line in f]

# define reference cycle
ref_cycle = '100'

# target start point and track length
starts = [[14.5,-167],
          [13.5,-164.5],
          [12.5,-162],
          [10.5,-160]]
track_length = 1000
direction = 'N'

# download and read ssha data for ascending tracks South of Hawaiian Ridge
Hawaii_SA_raw = Jason2.nested_dict()
for track in range(4):
   
    # define start and reference file:
    start = starts[track]
    reference_file = Jason2.download.find_track(start,ref_cycle,direction = direction)
    track_name = 'SA' + str(track + 1)
    
    # download, read data from each cycle
    for cycle in cycles:
        cycle_data = Jason2.download.read_track(start,track_length,ref_cycle,cycle,reference_file)
        if cycle_data == None:
            continue
        Hawaii_SA_raw[track_name][cycle] = cycle_data 

Copernicus

(NOTE: change directory names to where nc files are stored)

# set directories to retrieve nc files, initialize data structure
dir_names = ['D:/copernicus_current_data/{year}/{year}'.format(year=year) for year in np.arange(2008,2017)]
current_data = dict()

Hawaii_SA = Jason2.nested_dict()
for t_num,t_name in enumerate(Hawaii_SA_raw.keys()):
    
    # preprocess data 
    Hawaii_SA[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SA_raw[t_name],starts[t_num],distance_tol=15,masked_tol=5,supress_message=True)
    
    # find latitude and longitude range
    lats = np.array([[Hawaii_SA[t_name][cycle]['lat'][0],Hawaii_SA[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SA[t_name].keys()])
    lons = np.array([[Hawaii_SA[t_name][cycle]['lon'][0],Hawaii_SA[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SA[t_name].keys()])
    lat_range = [min(lats[:,0]), max(lats[:,1])]
    lon_range = [min(lons[:,0]), max(lons[:,1])]
    
    # download datasets from .nc directory (change directory name as necessary)
    time,lat,lon,u,v = Jason2.download.read_current_data(lat_range,lon_range,dir_names)
    current_data[t_name] = dict(zip(['time','lat','lon','u','v'],[time,lat,lon,u,v]))

ECCO version 4

(1) Use the bespoke Python package found here to read in $\frac{d \rho}{dr}$ datasets.

(2) Solve for buoyancy frequency using Jason2_py.find_N_from_DRHODR()

(3) Use Dedalus to solve Sturm-Liouville BVP for horizontal velocities of internal tides (see derivation). This can be done using script dedalus\solve_vertical_modes.py

For simplicity, we can just load the results from mode_plot.py from saved files:

Download datasets from saved files

# Jason-2
Hawaii_SA_raw = np.load('./test_data/Hawaii_SA_raw.npy',allow_pickle='TRUE').item()

# Copernicus
current_data = np.load('./test_data/current_dataSA.npy',allow_pickle='TRUE').item()

# ECCOv4
models = Jason2.nested_dict()
tracks = ['SA' + str(i+1) for i in range(4)]
for track in tracks:
    modes = 1 / np.genfromtxt('./test_data/eigenvectors/modes%s.csv'%track,delimiter = ',')[:,1]
    models[track]['M2 model']['wavenumber'] = modes[0]
    models[track]['S2 model']['wavenumber'] = modes[1]
    models[track]['eigenvectors'] = np.genfromtxt('./test_data\eigenvectors\%s.csv'%track,delimiter = ',')

Ascending Tracks

In this demonstration, we will analyze the 2D power spectra of principal semi-diurnal internal tides (M2 and S2) along tracks South of the Hawaiian Ridge (Figure 1).

Plot 2D power spectrum

Now that the data is downloaded, let's look compute the 2D power spectrum using Jason2_py.find_2D_power_spectrum():

# initialize plotting
Hawaii_SA = Jason2.nested_dict()

# starting coordinates
starts = [[14.5,-167],
          [13.5,-164.5],
          [12.5,-162],
          [10.5,-160]]


for t_num, t_name in enumerate(Hawaii_SA_raw.keys()):
    
    # preprocess dataset
    Hawaii_SA[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SA_raw[t_name],starts[t_num],distance_tol=15,masked_tol=5,supress_message=True)
    
    # find power spectrum
    fft_2D, dx, dt, nu_vector, f_vector = Jason2.find_2D_power_spectrum(Hawaii_SA[t_name])
    
    # find average lat and lon range
    lats = np.array([[Hawaii_SA[t_name][cycle]['lat'][0],Hawaii_SA[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SA[t_name].keys()])
    lons = np.array([[Hawaii_SA[t_name][cycle]['lon'][0],Hawaii_SA[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SA[t_name].keys()])
    lat_range = [np.mean(lats[:,0]), np.mean(lats[:,1])]
    lon_range = [np.mean(lons[:,0]), np.mean(lons[:,1])]
    
    # save info to models
    models[t_name]['data'].update({'fft 2D' : fft_2D,
                                  'dx': dx,
                                  'dt': dt,
                                  'wavenumbers' : nu_vector,
                                  'frequencies' : f_vector,
                                  'lat range' : lat_range,
                                  'lon range' : lon_range})

Below, we plot the surveyed Jason-2 tracks and their 2D power spectra:

# draw plots
track_dict = {'Track ' + str(i + 1) : pf.track_plt(models[track]['data'], 'Track ' + str(i + 1)) for i, track in enumerate(tracks)}
power_dict = {'Track ' + str(i + 1) : pf.fft_2D_plt(models[track]['data'], '') for i, track in enumerate(tracks)}
track_plot = hv.HoloMap(track_dict, kdims='Track').opts(infer_projection=True)
power_plot = hv.HoloMap(power_dict, kdims='Track')
(track_plot + power_plot).opts(title='Figure 1: 2D Power Spectra')

Models for Two Dimensional Fourier Transform

Anisotropic Wave

We begin model tides as sinusoids given by the function:

$$h(x,t) = A sin(2 \pi (k x + \omega t))$$

where $A$ is wave amplitude, $k$ is the tidal wavenumber and $\omega$ is the alias frequency. In our analysis, we minimize the effects of discontiuities by multiplying by a Hann window in both the space (x) and time (t). This window function is:

$$w(x,t) = \frac{1}{4}(1+cos(\frac{2 \pi x}{L})(1+cos(\frac{2 \pi t}{T}))rect(\frac{x}{L})rect(\frac{t}{T})$$

where $L$ is the length of the track and $T$ is the length of the time window. Finally, to account for the discretization of the signal, we multiply by a Dirac comb given by:

$$d(x,t) = \sum^\infty_{n=-\infty}\delta(t - n \Delta t)\sum^\infty_{n=-\infty}\delta(x - n \Delta x)$$

where $\Delta t$ and $\Delta x$ are the time and space intervals between sucessive measurements. In total, our signal function is $s(x,t) = h(x,t) \ w(x,t) \ d(x,t)$. We would like to generate an analytical function for the two dimensional Fourier transform, $\hat{s}(\nu,f) = \mathcal{F}[s(x,t)]$. Taking the fourier transforms of $h(x,t)$, $w(x,t)$, and $d(x,t)$, we have:

$$ \begin{align} \hat{h}(\nu,f) &= \frac{1}{2}A i \Big(\delta (\nu + k) \delta (f + \omega) - \delta (\nu - k) \delta (f - \omega) \Big) \\ \hat{w}(\nu,f) &= \frac{1}{4} \Big(\mathcal{F}[1 + cos(\frac{2 \pi x}{L})] * \mathcal{F}[rect(\frac{1}{L})]\Big)\Big(\mathcal{F}[1 + cos(\frac{2 \pi t}{T})] * \mathcal{F}[rect(\frac{t}{T})]\Big) \\ &= \frac{1}{4} \Big[ \Big(\delta(\nu) + \frac{1}{2} \delta (\nu - 1/L) + \frac{1}{2} \delta (\nu + 1/L)\Big) * L sinc(L \nu) \Big] \Big[\Big(\delta(f) + \frac{1}{2} \delta (f - 1/T) + \frac{1}{2} \delta (f + 1/T)\Big) * T sinc(T f) \Big] \\ &= \frac{1}{4}\Big[ L sinc(L \nu) + \frac{1}{2} L sinc(L \nu - 1) + \frac{1}{2} L sinc(L \nu + 1)\Big]\Big[ T sinc(T f) + \frac{1}{2} T sinc(T f - 1) + \frac{1}{2} T sinc(T f + 1)\Big] \\ \hat{d}(\nu,f) &= \frac{1}{\Delta x \Delta t} \sum^\infty_{n=-\infty}\delta(\nu - \frac{n}{\Delta x}) \sum^\infty_{n=-\infty}\delta(f - \frac{n}{\Delta t}) \end{align} $$

Thus, in total, we have the final analytical expression for the two-dimensional fourier transform of the signal:

$$ \begin{align} \hat{s}(\nu,f) &= \frac{A i}{8 \Delta x \Delta t} \\ & \times\Big[ \sum^\infty_{n=-\infty}\Big(L sinc(L(\nu - \frac{n}{\Delta x} + k)) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} + k) - 1) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} + k) + 1) \Big)\\ &\times \sum^\infty_{n=-\infty}\Big(T sinc(T(f - \frac{n}{\Delta t} + \omega)) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} + \omega) - 1) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} + \omega) + 1) \Big) \\ &- \sum^\infty_{n=-\infty}\Big(L sinc(L(\nu - \frac{n}{\Delta x} - k)) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} - k) - 1) + \frac{1}{2} L sinc(L(\nu - \frac{n}{\Delta x} - k) + 1) \Big)\\ &\times \sum^\infty_{n=-\infty}\Big(T sinc(T(f - \frac{n}{\Delta t} - \omega)) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} - \omega) - 1) + \frac{1}{2} T sinc(T(f - \frac{n}{\Delta t} - \omega) + 1) \Big) \Big] \\ \end{align} $$

The function Jason2_py.model.power_spectrum() computes the above result over a range of $n = [-50, 50]$. Let's use this to model the principal semi-diurnal tides, M2 (lunar) and S2 (solar).

# M2 and S2 tidal periods (hours)
periods = [12.4206012,12]
semidiurnals = ['M2 model', 'S2 model']
for t_num, t_name in enumerate(tracks):
    
    # unpack parameters from data
    track = models[t_name]['data']
    dx = track['dx']
    dt = track['dt']
    nu_vector = track['wavenumbers']
    f_vector = track['frequencies']
    L = dx * len(nu_vector)
    T = dt * len(f_vector)
    fft_2D = track['fft 2D']
    
    for i,tide in enumerate(semidiurnals):
        
        model = models[t_name][tide]
        k = model['wavenumber']
        
        # find alias frq from Jason-2 repeat period
        w = Jason2.helper.find_alias_frq(9.9156,periods[i]/24)
        model['alias frq'] = w
        
        # simulate spectrum, unit amplitude
        model_spectrum = Jason2.model.power_spectrum(1,L,nu_vector,k,T,w,f_vector,dx,dt)
    
        # fit model to predicted tidal peak, save values for plotting 
        model['fft 2D'] = Jason2.model.fit_models([[k,w]],nu_vector,f_vector,fft_2D, [model_spectrum]) * model_spectrum
        model['wavenumbers'] = nu_vector
        model['frequencies'] = f_vector
        model['frequencies'] = f_vector

Let's compare the model to the data (Figure 2):

np.seterr(divide='ignore')
# put plots into dict
_model_p = dict()
for i, track in enumerate(tracks):
    for tide in ['M2 model', 'S2 model']:
        _model_p[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track][tide],'model')

# make drop-down list
model_p = hv.HoloMap(_model_p, kdims=['Tide','Track'])

# plot overlay
(track_plot + 
 power_plot.opts(width =300, 
                 height = 300,
                 title='data',
                 colorbar=False) + 
 model_p.opts(width =300,
         height = 300,
         title='model')).opts(title='Figure 2')

We can also plot cross-sections centered at the estimated tidal peaks (Figure 3):

# initialize dict to hold plots
_f = dict()
_nu = dict()
for i, track in enumerate(tracks):
    
    # get data spectrum
    data = models[track]['data']
    for tide in ['M2 model', 'S2 model']:
        
        # get model spectrum
        model = models[track][tide]
        val_f = model['alias frq']
        val_nu = model['wavenumber']
        
        # make individual line plots
        mf = pf.line_plt(model,'f',val_nu).relabel('model')
        df = pf.line_plt(data,'f',val_nu).relabel('data')
        mnu = pf.line_plt(model,'nu',val_f).relabel('model')
        dnu = pf.line_plt(data,'nu',val_f).relabel('data')
        
        # make overlay dicts
        _f[(tide.replace(" model",""),'Track ' + str(i + 1))] = (df * mf).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
        _nu[(tide.replace(" model",""),'Track ' + str(i + 1))] = (dnu * mnu).opts(projection=ccrs.PlateCarree(),legend_position='top_left')

# # Make drop down list
f = hv.HoloMap(_f, kdims=['Tide','Track']).opts(infer_projection=True)
nu = hv.HoloMap(_nu, kdims=['Tide','Track']).opts(infer_projection=True)

# # plot layout
(f + nu).opts(title='Figure 3',shared_axes=False)

Anisotropic Wave with Frequency Doppler Shift

Comparing our M2 model to the data in frequency space, we find that the tidal peak is blue-shifted and broader than the model predicts. A possible explanation for this discrepancy is the effect of Doppler shifting by mean flows South of the Hawaiian Ridge. To examine this, we refine our model to include a time-dependent advection term using data from Copernicus. The Dopper shifted frequencies are calculated as:

$$\omega_t = \omega + {\overline{U}} \cdot \nu$$

where $\overline{U}$ is the current velocity averaged over the water column, $\nu$ is the tidal wavenumber, and $\omega$ is the alias frequency.

First, We can plot the surface velocities near the tracks (Figure 4):

_U = {'Track ' + str(i + 1) : pf.U_plt(current_data[track],models[track]['data'],"") for i, track in enumerate(tracks)}
U = hv.HoloMap(_U, kdims='Track').opts(infer_projection=True)
U.opts(title='Figure 4: Surface Current Magnitudes')

Thus, from the average surface velocity, we would expect the maximum magntiude of Doppler-shifting to be:

for i,track in enumerate(tracks):
    theta = np.arctan(np.diff(models[track]['data']['lat range']) / np.diff(models[track]['data']['lat range']))
    u = np.mean(current_data[track]['u'])
    v = np.mean(current_data[track]['v'])
    mag = np.sqrt(u**2 + v**2) * np.cos(theta)
    print(mag[0] * 24 * 3600 / 1000 * models[track]['M2 model']['wavenumber'], ' 1/day for track', str(i+1)) 

0.010002422659224665  1/day for track 1
0.01704400235270866  1/day for track 2
0.021981349014374017  1/day for track 3
0.03648486778777556  1/day for track 4

However, a better approximation would be to take the depth-averaged velocities, which can be found by solving for the vertical modes, as derived below.

Calculation of $k$ and $ U$

The horizontal velocities of internal tides for a rigid-lid, flat-bottom ocean are described by the following Sturm-Liouville boundary value problem (1).:

$$\frac{d}{dz}\Big(\frac{1}{N^2}\frac{dF}{dz}\Big) + \frac{1}{c^2}F = 0, \\ \frac{dF}{dz} = 0 \ \ \ \ \text{at} \ \ \ \ z = 0 \ \ \ \ \text{ and} \ \ \ \ z= -H$$

where $H$ is ocean depth and $N$ is the buoyancy frequency. The eigenvalue solutions (vertical modes), $F_n$, have associated eigenvalues of $-1/c_n^2$, and the tidal wavneumber is calculated using the dispersion relationship:

$$k_n = \frac{(\omega^2 - f^2)^{1/2}}{c_n}$$

where $f$ is the coriolis parameter and $\omega$ is the tidal frequency (e.g. reported here). $F_n$ are solved using the dedalus framework with stratification from ECCO version 4. Estimates of current velocity are done by assuming horizontal velocities are given by

$$U = A_1 F_1 + A_2 F_2, \\ U = 0 \ \ \ \ \text{at} \ \ \ \ z = -H \ \ \ \ \text{and} \ \ \ \ U = U_1 \ \ \ \ \text{at} \ \ \ \ z = 0$$

where $U_1$ is the current velocity reported by Copernicus projected on the Jason-2 track. In our model, Jason2_py.model.doppler_power_spectrum(), we simulate a power spectrum for each Jason-2 pass using the $\omega_t$ from that day and report the average of all simulations.

for t_num, t_name in enumerate(tracks):
    
    data = models[t_name]['data']
    eigenvector = models[t_name]['eigenvectors'][2,:]
    c_data = current_data[t_name]
    for tide in ['M2 model', 'S2 model']:
        print('simulating track ' + ''.join([str(t_num + 1), ' ', tide]))
        model = models[t_name][tide]
    
        # evaluate doppler model
        args = (data['lat range'],
               data['lon range'],
               data['frequencies'],
               data['wavenumbers'],
               model['alias frq'],
               data['dx'] * len(data['wavenumbers']),
               model['wavenumber'],
               data['dt'] * len(data['frequencies']),
               data['dx'],
               data['dt'],
               data['fft 2D'])
        
        w_t,power_spec  = Jason2.model.doppler_power_spectrum(eigenvector,c_data,Hawaii_SA[t_name],args)
        
        # fit model to tidal wavenumber and alias frequency
        coeff = Jason2.model.fit_models([[args[6],np.mean(w_t)]],args[3],args[2],args[-1], [power_spec])
        fft_2D = coeff * power_spec

        # store values
        models[t_name]['doppler ' + tide] = {'wavenumber' : model['wavenumber'],
                                             'alias frq'  : model['alias frq'],
                                             'shifted frequencies' : w_t,
                                             'fft 2D' : fft_2D,
                                             'wavenumbers' : data['wavenumbers'],
                                             'frequencies' : data['frequencies']}

simulating track 1 M2 model
simulating track 1 S2 model
simulating track 2 M2 model
simulating track 2 S2 model
simulating track 3 M2 model
simulating track 3 S2 model
simulating track 4 M2 model
simulating track 4 S2 model

Calculating the $ \omega_t$ over all Jason-2 cycles, we have the following distribution (Figure 5):

_wt = dict()
for i, track in enumerate(tracks):
    frequencies, edges = np.histogram(models[track]['doppler M2 model']['shifted frequencies'], 20)
    title = 'Alias Frequency: %.5f 1/day'%models[track]['M2 model']['alias frq']
    _wt['Track ' + str(i + 1)] = hv.Histogram((edges, frequencies)).opts(xlabel = '𝜔ₜ (1/day)',
                                                                         title = title,
                                                                         fontsize = {'labels':12})
wt = hv.HoloMap(_wt, kdims='Track').opts(title='Figure 5')
wt

We can plot the power spectra of the doppler simulations and their cross-sections (Figure 6 and 7):

# put plots into dict
_doppler = dict()
for i, track in enumerate(tracks):
    for tide in ['M2 model', 'S2 model']:
        _doppler[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track]['doppler ' + tide],'doppler model')

# make drop-down list
doppler = hv.HoloMap(_doppler, kdims=['Tide','Track']).opts(width =300,height = 300)

# plot overlay
p1 = power_plot.opts(width = 250, height = 300,title='data', colorbar=False)
p2 = model_p.opts(width = 250,height = 300,title='model',colorbar=False)
(p1 + p2 + doppler).opts(title='Figure 6')

# initialize dict to hold plots
_frq = dict()
_om = dict()
for i, track in enumerate(tracks):
    
    # get data spectrum
    data = models[track]['data']
    
    for tide in ['M2 model', 'S2 model']:
        
        # get model spectrum
        model = models[track][tide]
        dmodel = models[track]['doppler ' + tide]
        val_f = model['alias frq']
        val_nu = model['wavenumber']
        
        # make individual line plots
        df = pf.line_plt(data,'f',val_nu).relabel('data')
        mf = pf.line_plt(model,'f',val_nu).relabel('model')
        dmf = pf.line_plt(dmodel,'f',val_nu).relabel('doppler model')
        dnu = pf.line_plt(data,'nu',val_f).relabel('data')
        mnu = pf.line_plt(model,'nu',val_f).relabel('model')
        dmnu = pf.line_plt(dmodel,'nu',val_f).relabel('doppler model')
        
        # make overlay dicts
        _frq[(tide.replace(' model',''),'Track ' + str(i + 1))] = (df * mf * dmf).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
        _om[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dnu * mnu * dmnu).opts(projection=ccrs.PlateCarree(),legend_position='top_left')

# Make drop down list
frq = hv.HoloMap(_frq, kdims=['Tide','Track']).opts(infer_projection=True)
om = hv.HoloMap(_om, kdims=['Tide','Track']).opts(infer_projection=True)

# plot layout
(frq.opts(width=450) + om.opts(width=400)).opts(title='Figure 7',shared_axes=False)

Below, we have a summary of mean square error:

# ranges (track x frequency range x wavenumber range)
ranges = np.array([[[.008,.03],[0.0055,0.008]],
                  [[.01,.03],[.003,.0085]],
                  [[0.000,.04,],[.005,.0075]],
                  [[.01,.025],[.005,.009]]])

# calculate mean square error
mse = np.zeros((4,2,2))
_range = dict()
for t_num,t_name in enumerate(tracks):
    
    # get ranges
    f_range = ranges[t_num,0,:].tolist()
    nu_range = ranges[t_num,1,:].tolist()
    _range['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'],'').opts(xlim=tuple(nu_range),
                                                                                     ylim=tuple(f_range))
    # get data spectrum
    power_spectrum = models[t_name]['data']
    f_vector = power_spectrum['frequencies']
    nu_vector = power_spectrum['wavenumbers']
    for tide_num, tide in enumerate(['M2 model', 'S2 model']):
        
        # get model specta
        model = models[t_name][tide]['fft 2D']
        dmodel = models[t_name]['doppler ' + tide]['fft 2D']
        
        # solve for mean square error
        mse[t_num,tide_num,0] = Jason2.find_square_error(f_range, 
                                                         nu_range,
                                                         f_vector,
                                                         nu_vector,
                                                         power_spectrum['fft 2D'],
                                                         model)
        mse[t_num,tide_num,1]= Jason2.find_square_error(f_range,
                                                        nu_range,
                                                        f_vector,
                                                        nu_vector,
                                                        power_spectrum['fft 2D'],
                                                        dmodel)
# output
range_plt = hv.HoloMap(_range, kdims='Track').opts(framewise=True,
                                                  width=400,height=300)
d = dict()
d['M2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
                                data=mse[:,0,:],
                                index = ['Track %.f'%(i+1) for i in range(4)])

d['S2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
                                data=mse[:,1,:],
                                index = ['Track %.f'%(i+1) for i in range(4)])

pd.set_option('display.float_format', '{:.2E}'.format)
df = pd.concat(d, axis=1,names=['Tidal Constituent:'])

print('Table 1: Mean Square Error')
display(HTML(df.to_html()))
range_plt.opts(title='Figure 8: Semi-Diurnal Peaks')

Table 1: Mean Square Error

Tidal Constituent:	M2		S2
	Model	Doppler Model	Model	Doppler Model
Track 1	2.02E-04	8.76E-05	1.30E-04	1.25E-04
Track 2	5.11E-04	4.92E-04	5.23E-04	5.01E-04
Track 3	1.29E-04	3.51E-05	1.24E-04	5.28E-05
Track 4	2.39E-05	2.19E-05	2.93E-05	2.98E-05

Thus over the ranges selected, we see that the Doppler Model reduces the mean square error for all but track 4. Let's now see if we can pick up the same tidal signals from the descending tracks.

Descending Tracks

Click here to load data (or run next two code cells).

Load Jason-2 Data

# get existing cycles:
f = open("cycles.txt", "r")
cycles = [line.strip()[1:4] for line in f]

# define reference cycle
ref_cycle = '100'

# target start point and track length
starts = [[22.450964, -168.03407],
          [21.74352, -164.891574],
          [20.98662, -161.731304],
          [18.619012, -160.738939]]
track_length = 1000
direction = 'S'

# download and read ssha data for ascending tracks South of Hawaiian Ridge
Hawaii_SD_raw = Jason2.nested_dict()
for track in range(4):
   
    # define start and reference file:
    start = starts[track]
    reference_file = Jason2.download.find_track(start,ref_cycle,direction = direction)
    track_name = 'SD' + str(track + 1)
    
    # download, read data from each cycle
    for cycle in cycles:
        cycle_data = Jason2.download.read_track(start,track_length,ref_cycle,cycle,reference_file)
        if cycle_data == None:
            continue
        Hawaii_SD_raw[track_name][cycle] = cycle_data 

Load Copernicus Data (set directory name as necessary)

# set directories to retrieve nc files, initialize data structure
dir_names = ['D:/copernicus_current_data/{year}/{year}'.format(year=year) for year in np.arange(2008,2017)]
currents_SD = dict()

Hawaii_SD = Jason2.nested_dict()
lat_range = np.zeros((4,2))
lon_range = np.zeros((4,2))
for t_num, t_name in enumerate(Hawaii_SD_raw.keys()):
    
    # preprocess data
    Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)
    
    # find latitude and longitude range
    lats = np.array([[Hawaii_SD[t_name][cycle]['lat'][0],Hawaii_SD[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SD[t_name].keys()])
    lons = np.array([[Hawaii_SD[t_name][cycle]['lon'][0],Hawaii_SD[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SD[t_name].keys()])
    lat_range = [min(lats[:,1]), max(lats[:,0])]
    lon_range = [min(lons[:,0]), max(lons[:,1])]
    
    # download datasets from .nc directory (change directory name as necessary)
    time,lat,lon,u,v = Jason2.download.read_current_data(lat_range,lon_range,dir_names)
    currents_SD[t_name] = dict(zip(['time','lat','lon','u','v'],[time,lat,lon,u,v]))

Load Descending Track Data from saved files

# Jason-2
Hawaii_SD_raw = np.load('./test_data/Hawaii_SD_raw.npy',allow_pickle='TRUE').item()

# Copernicus
current_data = np.load('./test_data/current_dataSD.npy',allow_pickle='TRUE').item()

# ECCOv4
tracksd = ['SD' + str(i+1) for i in range(4)]
for track in tracksd:
    modes = 1 / np.genfromtxt('./test_data\eigenvectors\modes%s.csv'%track,delimiter = ',')[:,1]
    models[track]['M2 model']['wavenumber'] = modes[0]
    models[track]['S2 model']['wavenumber'] = modes[1]
    models[track]['eigenvectors'] = np.genfromtxt('./test_data\eigenvectors\%s.csv'%track,delimiter = ',')  

Angular Dependence of Tidal Wavenumber

Unlike the ascending tracks, we assume that waves do not propagate collinearly with the Jason-2 tracks. Thus, the effective wavenumbers used to describe the observed 2D power spectra are given as:

$$k_{eff} = k \ cos(\theta)$$

where $k$ is the expected tidal wavenumber, and $\theta$ is the angle between the tide beam and the track. Since we wnat to model the strong semi-dirunal signals found in the ascending tracks, we have chosen 4 descending tracks whose midpoints are closest to the respective midpoints of the 4 ascending tracks. This choice was made because the Hann window places greater weight on values near the center of each track.

Below, we calculate $\theta$ for each track by simply taking the angle between each ascending-descending track pair:

Hawaii_SD = Jason2.nested_dict()
_trkd = dict()
_trk = dict()
tracksd = ['SD' + str(i+1) for i in range(4)]
starts = [[22.450964, -168.03407],
          [21.74352, -164.891574],
          [20.98662, -161.731304],
          [18.619012, -160.738939]]

thetas = [0] * 4
for t_num, t_name in enumerate(tracksd):
    
    # preprocess data
    Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)  
    
    # find average lat and lon range
    lats = np.array([[Hawaii_SD[t_name][cycle]['lat'][0],Hawaii_SD[t_name][cycle]['lat'][-1]] for cycle in Hawaii_SD[t_name].keys()])
    lons = np.array([[Hawaii_SD[t_name][cycle]['lon'][0],Hawaii_SD[t_name][cycle]['lon'][-1]] for cycle in Hawaii_SD[t_name].keys()])
    lat_range = [np.mean(lats[:,0]), np.mean(lats[:,1])]
    lon_range = [np.mean(lons[:,0]), np.mean(lons[:,1])]
    
    # save info to models
    models[t_name]['data'].update({'lat range' : lat_range,
                                  'lon range' : lon_range})
    
    # find track angles:
    data_a = models['SA' + str(t_num + 1)]['data']
    data_d = models[t_name]['data']
    v_a = [np.diff(data_a['lat range'])[0],np.diff(data_a['lon range'])[0]]
    v_d = [np.diff(data_d['lat range'])[0],np.diff(data_d['lon range'])[0]]
    thetas[t_num] = np.arccos(np.dot(v_d/np.linalg.norm(v_d),v_a/np.linalg.norm(v_a)))
    
    # create plot dicts
    _trkd['Track ' + str(t_num+1)] = pf.track_plt(models[t_name]['data'],"Track %.f"%(t_num+1))
    _trk['Track ' + str(t_num+1)] = (track_dict['Track ' + str(t_num+1)] * _trkd['Track ' + str(t_num+1)]).opts(projection = ccrs.PlateCarree(),
                                                                                                                title='Figure 9: θ = %.3f'%thetas[t_num]) 

# make drop down plots
trk = hv.HoloMap(_trk,kdims=['Track'])
trk

Plotting the 2D power spectra of the descending tracks:

_powerd = dict()
for t_num, t_name in enumerate(tracksd):
    
    # preprocess data
    Hawaii_SD[t_name] = Jason2.preprocess.preprocess_dataset(Hawaii_SD_raw[t_name],starts[t_num],masked_tol=5)    
    
    # find power spectrum
    fft_2D, dx, dt, nu_vector, f_vector = Jason2.find_2D_power_spectrum(Hawaii_SD[t_name])
    
    # save info to models
    models[t_name]['data'].update({'fft 2D' : fft_2D,
                                  'dx': dx,
                                  'dt': dt,
                                  'wavenumbers' : nu_vector,
                                  'frequencies' : f_vector})
    
    # Make heatmaps
    _powerd['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'], '')
    

# make drop-down plots
powerd = hv.HoloMap(_powerd, kdims=['Track'])
(trk.opts(title='',infer_projection=True,height=200,width=200) +
 power_plot.opts(title='Ascending',height=300,width=300) +
 powerd.opts(title='Descending',height=300,width=300)).opts(title='Figure 10: 2D Power Spectra')

As expected, we are fitting to the peaks where the tidal wavenumber and alias frequency are opposite in sign. Now, let's turn to modeling. Plotting the surface currents we see similar magnitudes as for the ascending tracks, so we expect the Doppler shifting of the alias frequency to be a factor.

_Ud = {'Track ' + str(i + 1) : pf.U_plt(current_data[track],models[track]['data'],"") for i, track in enumerate(tracksd)}
Ud = hv.HoloMap(_Ud, kdims='Track').opts(infer_projection=True)
Ud.opts(title='Figure 11: Mean Flow at Surface')

Let's compute both the simple and Doppler model, and update our dictionaries:

# M2 and S2 tidal periods (hours)
periods = [12.4206012,12]
semidiurnals = ['M2 model', 'S2 model']

for t_num, t_name in enumerate(tracksd):
    
    # get data
    c_data = current_data[t_name]
    data = models[t_name]['data']
    eigenvector = models[t_name]['eigenvectors'][2,:] # 2nd vertical mode
    dx = data['dx']
    dt = data['dt']
    nu_vector = data['wavenumbers']
    f_vector = data['frequencies']
    L = dx * len(nu_vector)
    T = dt * len(f_vector)
    fft_2D = data['fft 2D']
    
    for i,sdiurnal in enumerate(semidiurnals):
        
        model = models[t_name][sdiurnal]
        k = model['wavenumber']
        
        # multiply k by track angle
        k_eff = k * np.cos(thetas[t_num])
        
        # find alias frq from Jason-2 repeat period
        w = Jason2.helper.find_alias_frq(9.9156,periods[i]/24)
        model['alias frq'] = w
        
        
        #-------Model 1: No Doppler shift--------
        print('simulating track ' + ''.join([str(t_num + 1), ' ', sdiurnal]))
        model_spectrum = Jason2.model.power_spectrum(1,L,nu_vector,k_eff,T,w,f_vector,dx,dt)
        model['fft 2D'] = Jason2.model.fit_models([[k_eff,w]],nu_vector,f_vector,fft_2D, [model_spectrum]) * model_spectrum
        model['wavenumbers'] = nu_vector
        model['frequencies'] = f_vector
        model['frequencies'] = f_vector
        model['effective wavenumber'] = k_eff
        
        #-------Model 2: Doppler Shift---------
        args = (data['lat range'],
               data['lon range'],
               data['frequencies'],
               data['wavenumbers'],
               model['alias frq'],
               data['dx'] * len(data['wavenumbers']),
               k_eff,
               data['dt'] * len(data['frequencies']),
               data['dx'],
               data['dt'],
               data['fft 2D'])
        
        w_t,power_spec  = Jason2.model.doppler_power_spectrum(eigenvector,c_data,Hawaii_SD[t_name],args)
        
        # fit model to tidal wavenumber and average Doppler-shifted frequency
        coeff = Jason2.model.fit_models([[args[6],np.mean(w_t)]],args[3],args[2],args[-1], [power_spec])
        fft_2D = coeff * power_spec
        models[t_name]['doppler ' + sdiurnal] = {'wavenumber' : model['wavenumber'],
                                                'effective wavenumber' : k_eff,
                                                'alias frq'  : model['alias frq'],
                                                'shifted frequencies' : w_t,
                                                'fft 2D' : fft_2D,
                                                'wavenumbers' : data['wavenumbers'],
                                                'frequencies' : data['frequencies']}

simulating track 1 M2 model
simulating track 1 S2 model
simulating track 2 M2 model
simulating track 2 S2 model
simulating track 3 M2 model
simulating track 3 S2 model
simulating track 4 M2 model
simulating track 4 S2 model

# put plots into dict
_modeld = dict()
_dmodeld = dict()
for i, track in enumerate(tracksd):
    for tide in ['M2 model', 'S2 model']:
        _modeld[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track][tide],'model')
        _dmodeld[(tide.replace(" model",""),'Track ' + str(i + 1))] = pf.fft_2D_plt(models[track]['doppler ' + tide],'model')
# make drop-down list
modeld = hv.HoloMap(_modeld, kdims=['Tide','Track'])
dmodeld = hv.HoloMap(_dmodeld, kdims=['Tide','Track'])
trkd = hv.HoloMap(_trkd, kdims =['Track']).opts(infer_projection=True)

# plot overlay
(powerd.opts(width =250, 
                 height = 300,
                 title='data',
                 colorbar=False) + 
 modeld.opts(width =250,
         height = 300,
         title='model',
         colorbar=False) + 
 dmodeld.opts(width =300,
         height = 300,
         title='Doppler model')).opts(title='Figure 12: 2D models')

(powerd + modeld).opts(shared_axes=True)

# initialize dict to hold plots
_frqd = dict()
_omd = dict()
for i, track in enumerate(tracksd):
    
    # get data spectrum
    data = models[track]['data']
    
    for tide in ['M2 model', 'S2 model']:
        
        # get model spectrum
        model = models[track][tide]
        dmodel = models[track]['doppler ' + tide]
        val_f = model['alias frq']
        val_nu = model['effective wavenumber']
        
        # make individual line plots
        dfd = pf.line_plt(data,'f',val_nu).relabel('data')
        mfd = pf.line_plt(model,'f',val_nu).relabel('model')
        dmfd = pf.line_plt(dmodel,'f',val_nu).relabel('doppler model')
        dnud = pf.line_plt(data,'nu',val_f).relabel('data')
        mnud = pf.line_plt(model,'nu',val_f).relabel('model')
        dmnud = pf.line_plt(dmodel,'nu',val_f).relabel('doppler model')
        
        # make overlay dicts
        _frqd[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dfd * mfd * dmfd).opts(projection=ccrs.PlateCarree(),legend_position='top_left')
        _omd[(tide.replace(' model',''),'Track ' + str(i + 1))] = (dnud * mnud * dmnud).opts(projection=ccrs.PlateCarree(),legend_position='top_left')

# Make drop down list
frqd = hv.HoloMap(_frqd, kdims=['Tide','Track']).opts(infer_projection=True)
omd = hv.HoloMap(_omd, kdims=['Tide','Track']).opts(infer_projection=True)

# plot layout
(frqd.opts(width=450) + omd.opts(width=400)).opts(title='Figure 13',shared_axes=False)

Descending Track Models: Mean Square Error

# ranges (track x frequency range x wavenumber range)
rangesd = np.array([[[.01,.03],[-.008,-.0035]],
                  [[.01,.03],[-.008,-.0035]],
                  [[0.01,.03,],[-.006,-.003]],
                  [[.015,.025],[-.008,-.004]]])

# calculate mean square error
msed = np.zeros((4,2,2))
_ranged = dict()
for t_num,t_name in enumerate(tracksd):
    
    # get ranges
    f_range = rangesd[t_num,0,:].tolist()
    nu_range = rangesd[t_num,1,:].tolist()
    _ranged['Track ' + str(t_num + 1)] = pf.fft_2D_plt(models[t_name]['data'],'').opts(xlim=tuple(nu_range),
                                                                                      ylim=tuple(f_range))
    # get data spectrum
    power_spectrum = models[t_name]['data']
    f_vector = power_spectrum['frequencies']
    nu_vector = power_spectrum['wavenumbers']
    for tide_num, tide in enumerate(['M2 model', 'S2 model']):
        
        # get model specta
        model = models[t_name][tide]['fft 2D']
        dmodel = models[t_name]['doppler ' + tide]['fft 2D']
        
        # solve for mean square error
        msed[t_num,tide_num,0] = Jason2.find_square_error(f_range, 
                                                         nu_range,
                                                         f_vector,
                                                         nu_vector,
                                                         power_spectrum['fft 2D'],
                                                         model)
        msed[t_num,tide_num,1]= Jason2.find_square_error(f_range,
                                                        nu_range,
                                                        f_vector,
                                                        nu_vector,
                                                        power_spectrum['fft 2D'],
                                                        dmodel)
# output
range_pltd = hv.HoloMap(_ranged, kdims='Track').opts(framewise=True,
                                                  width=400,height=300)
dd = dict()
dd['M2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
                                data=msed[:,0,:],
                                index = ['Track %.f'%(i+1) for i in range(4)])

dd['S2'] = pd.DataFrame(columns=['Model', 'Doppler Model',],
                                data=msed[:,1,:],
                                index = ['Track %.f'%(i+1) for i in range(4)])

pd.set_option('display.float_format', '{:.2E}'.format)
dfd = pd.concat(dd, axis=1,names=['Tidal Constituent:'])

print('Table 2: Mean Square Error')
display(HTML(dfd.to_html()))
range_pltd.opts(title='Figure 8: Semi-Diurnal Peaks')

Table 2: Mean Square Error

Tidal Constituent:	M2		S2
	Model	Doppler Model	Model	Doppler Model
Track 1	2.69E-04	1.13E-04	2.07E-04	9.04E-05
Track 2	1.01E-03	7.84E-04	8.94E-04	6.67E-04
Track 3	1.42E-04	6.52E-05	1.15E-04	7.35E-05
Track 4	3.12E-05	2.12E-05	3.01E-05	2.50E-05

References

(1) https://oceanphysics.files.wordpress.com/2019/08/jpo-d-18-0272.1-1.pdf